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XMDS2 is a cross-platform, GPL-licensed, open source package for numerically integrating initial value problems that range from 
a single ordinary differential equation up to systems of coupled stochastic partial differential equations. The equations are described 
in a high-level XML-based script, and the package generates low-level optionally parallelised C++ code for the efficient solution 
of those equations. It combines the advantages of high-level simulations, namely fast and low-error development, with the speed, 
portability and scalability of hand-written code. XMDS2 is a complete redesign of the XMDS package, and features support for a 
much wider problem space while also producing faster code. 
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i 1 Program summary 

Manuscript Title: XMDS2: Fast, scalable simulation of coupled 
I stochastic partial differential equations 
Q" 1 Authors: Graham R. Dennis, Joseph J. Hope, Mattias T. Johnsson 
£h Program Title: XMDS 2 
Journal Reference: 
Catalogue identifier: 

Licensing provisions: GNU Public License 2 
Programming language: Python and C++. 
^^Computer: Any computer with a Unix-like system, a C++ compiler 
f-*j and Python 

| Operating system: Any Unix-like system; developed under Mac OS X 
l— ~~ 1 and GNU/Linux 

RAM: Problem dependent (roughly 50 bytes per grid point) 
Number of processors used: Up to the minimum of the number of 
points in each of the first two dimensions 

Keywords: Initial value problems, differential equations, stochastic 
partial differential equations 

Classification: 4.3 Differential Equations, 6.5 Software including Par- 
allel Algorithms 

External routines/libraries: The external libraries required are 
problem-dependent. Uses FFTW3 Fourier transforms (used only 
^ for FFT-based spectral methods), dSFMT random number generation 
. ,—i (used only for stochastic problems), MPI message-passing interface 
(used only for distributed problems), HDF5, GNU Scientific Library 
(used only for Bessel-based spectral methods) and a BLAS implemen- 
tation (used only for non-FFT-based spectral methods). 
Nature of problem: General coupled initial-value stochastic partial dif- 
ferential equations 

Solution method: Spectral method with method-of-lines integration 
Running time: Determined by the size of the problem 
Web site: http : / / www . xmds . org 
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1. Introduction 

The integration of a system of variables from a set of initial 
conditions is one of the most widely performed tasks in quanti- 
tative simulation. Numerical integration is typically performed 
in one of two different styles: high-level methods using general 
software tools, or low-level methods using bespoke hand-tuned 
source code. The high-level approach requires much less code, 
and is therefore fast to develop and comparatively free of cod- 
ing errors. However, the low-level approach can provide dra- 
matic and necessary performance improvements, can utilise the 
full capacity of the computing platform for which it is devel- 
oped, and is more customisable. XMDS2 is a software package 
whose aim is to provide the key benefits of both approaches (TJ. 

The purpose of XMDS2 is to simplify the process of creat- 
ing simulations that solve systems of initial-value partial and 
ordinary differential equations. Instead of going through the 
error-prone process of hand-writing thousands of lines of code, 
XMDS2 enables problems to be described in a simple XML 
format. From this XML description XMDS2 generates a C++ 
simulation that solves the problem using fast algorithms. The 
code generated by XMDS2 is typically as fast as, or faster than, 
hand-written code, but by using XMDS2 the time taken to pro- 
duce the simulation is significantly reduced. 

XMDS2 can be used to simulate almost any set of (coupled) 
(partial) (stochastic) differential equations in any number of di- 
mensions. It can input and output data in a range of data for- 
mats, produce programs that can take command-line arguments, 
and produce parallelised code suitable for either modern com- 
puter architectures or distributed clusters. 

Aside from innumerable low-level libraries and high-level 
packages for numerical integration, there have also been mul- 
tiple previous attempts to automate or semi-automate the pro- 
cess of coding low-level numerical simulations, such as ||2][3). 
Rather than provide 'shell code' that can be edited, an XMDS2 
script is effectively a self-contained language written in XML 
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which is used to generate a fast C++ simulation. 

The first version of XMDS was released in 1997 as an open- 
source software package written in C++ which could simulate 
a class of stochastic partial differential equations [4]. Over the 
next decade, it was then expanded in scope and features by a 
growing group of developers. While most of these developers 
came from the fields of quantum optics and atom optics, where 
its ability to integrate stochastic equations is particularly per- 
tinent, the package slowly gained wider popularity. In 2008, 
the decision was taken to completely rewrite the package in 
Python (although still generating low-level C++ code), with a 
re-engineered structure that would allow it to address a much 
wider problem space. XMDS2 was released in 2010, and has 
recently received its first major update, along with extensive 
documentation, installers, and an examples library. 

Citing only a few examples, XMDS and increasingly 
XMDS2 have been used in the fields of quantum atom optics 
15] El, quantum optics QUI, quantum control [9|, predator- 
prey dynamics IflOll and ecology IfTTl [12l [T3l . 

2. Problem class 

XMDS2 solves systems of initial-value differential equa- 
tions. Each differential equation can: 

1 . have an arbitrary number of dimensions, which may differ 
from that of other differential equations in the system, 

2. involve integrals of quantities in the differential system, or 

3. include stochastic elements either in initial conditions and 
filters, or in the dynamical equations themselves. 

As an example, property 1 means that XMDS2 can solve sys- 
tems in which a partial differential equation is coupled to an or- 
dinary differential equation. Property 2 allows the evolution of 
the ordinary differential equation to depend upon moments of 
the partial differential equation. Property 2 also allows partial 
differential equations to depend non-locally on system quanti- 
ties. Property 3 permits the integration of systems of stochastic 
(partial) differential equations, which are typically written us- 
ing either Gaussian noise (via a Wiener process), or a Poisso- 
nian noise in which the system changes state in a discontinuous 
way. 

XMDS2 uses spectral methods lfl4l . which induce two re- 
strictions on the problem space. The first is that the geometry 
of the simulation domain must be a tensor product of lattices 
in each dimension (see Figure[T]l. The second restriction is that 
the boundary conditions must be compatible with the spectral 
method used. XMDS2 currently supports periodic, even, odd 
and zero boundary conditions. Spectral methods allow the ap- 
proximation of spatial derivatives with 'exponential' accuracy 
In addition, the restriction to tensor product 



(see Section 3.1 



lattices affords significant computational savings which will be 
discussed later. The disadvantage is that XMDS2 cannot be 
used to solve problems on arbitrary-shaped domains as is pos- 
sible using finite-element methods. This is not a significant lim- 
itation for a wide class of problems, as the system is often con- 
strained to evolve within a finite domain. Quantum atom optics 
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Figure 1: An example of a tensor product lattice. The specific unequal lat- 
tice spacing in the x direction is due to the fact that the basis functions in that 
dimension have been chosen to be Hermite-Gauss. 



problems, for example, have trapping terms in the differential 
equation that cause the solution to be non-negligible over a fi- 
nite domain. 

The use of spectral methods means that XMDS2 represents 
the solution as a linear combination of global basis functions 
that extend over the entire domain. This is an accurate represen- 
tation for solutions which are smooth. Problems which contain 
shocks or other spatial discontinuities (including discontinuous 
derivatives) are better served by local methods such as finite 
difference or finite element methods. 

Subject to these caveats, XMDS2 is applicable to a broad 
problem class and employs efficient and accurate algorithms for 
the solution of these problems. 



3. Algorithms employed 

XMDS2 employs efficient algorithms in its generated simu- 
lations. These include: 



1. 

2. 

3. 
4. 
5. 



Spectral methods for computing spatial derivatives, 

Fast spatial-to-spectral transforms including FFTs and 

parity-exploiting matrix transforms, 

Distributed memory parallelism, 

Gaussian quadrature for spatial integration, 

Method-of-lines explicit temporal integration schemes, 

and 

Interaction picture methods for exactly solving linear parts 
of the problem. 



3.1. 



The spectral method 

XMDS2 spatially discretises the problem by applying the 
spectral method [ 14] . This method decomposes the solution as 
a weighted sum of a finite set of orthonormal basis functions. 
For example, the quantity f(x,y) is represented as 



f(x,y) = 2^ F, hm X n (x)Y m (y), 



(1) 
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where F n m is a matrix of coefficients, X„(x) is the n th basis 
function for the x dimension, and Y,„(y) is the m th basis func- 
tion for the y dimension. The coefficients F„ m fully describe 
the solution and are the spectral representation of the solution. 
Typically in XMDS2, the number of basis functions is equal to 
the number of grid points in each dimension. In this case, the 
spectral representation is equivalent to the spatial representa- 
tion f(xi,yj), the values of the solution at the grid points. The 
two representations are linked by the linear transformation ([T} 
and its inverse. 

Spectral methods approximate spatial derivatives using the 
decomposition ([TJ and using analytic expressions for the deriva- 
tives of the basis functions, 



d p d q 
dxP dy 



f(x,y) = Y j F n , 



dPX„(x) dPYJy) 



dxP 



dyi 



(2) 



Spatial derivatives approximated in this manner are 'expo- 
nentially' accurate. In general, an optimal M-point method 
to calculate a £-order derivative of a function will have error 
0{h M ~ k ) where h is the grid-point spacing. As h oc 1 /N where 
N is the number of grid points, such a method will converge 
like 0(l/N M ~ k ) for a &-order derivative. In spectral methods 
the value of the solution at all grid points is used when com- 
puting spatial derivatives, hence M — N. In this case M in- 
creases as the number of grid points N increases resulting in a 
method whose order effectively increases as the number of grid 
points increases. The asymptotic error of a spectral method is 
0(l/N N ~ k ) which converges exponentially. 

The basis functions are typically chosen to make part of the 
differential equation diagonal in the spectral basis. XMDS2 
supports the following spectral methods for each dimension: 

• Fourier modes (complex exponentials), 

X„(x) = e ik " x . 

This method imposes periodic boundary conditions. The 
basis functions are eigenfunctions of the Cartesian spatial 
derivative operator. This is a general purpose method. 

• Cosine / sine functions, 

X„(x) = cos(k„x) or X„ = sin(fc„x). 

These methods impose even and odd boundary conditions 
respectively at the ends of the domain. The basis functions 
are eigenfunctions of the Laplacian in Cartesian coordi- 
nates. This method is useful when the problem has even 
or odd reflection symmetry about a plane. 

• 'Cylindrical' Bessel functions, 

R n (r) = J,„(k n r), 

where J,„(r) is the order-m Bessel function of the first kind. 
This method imposes analytic boundary conditions at the 
origin and zero Dirichlet boundary conditions at the outer 
boundary. The basis functions are eigenfunctions of the 
radial component of the Laplacian in cylindrical coordi- 
nates. This method is useful for problems with rotational 
symmetry. See 1151 for more details. 



'Spherical' Bessel functions, 



R n (r) 



2r J l+ i(k n r). 



This method imposes analytic boundary conditions at the 
origin and zero Dirichlet boundary conditions at the outer 
boundary. The basis functions are eigenfunctions of the ra- 
dial component of the Laplacian in spherical coordinates. 
This method is useful for problems with spherical symme- 
try. 

Hermite-Gauss functions, 

iff n (x) = {2 n n\o-y[ny m e- xLl1 ^ ' H„((rx), where: 
d" 



dx n 

This method requires that the solution decay as e~ x ^ 2tT in 
the limit x — > ±oo. The basis functions are eigenfunctions 
of the Schrodinger equation for the harmonic oscillator: 

H 2 d 2 if/„ 1 / 1 \ 

~2m~dx^ + 2 mc ° 2x2 ^" ( - x " > = haJ \ n + 2 j^"^' ^ 



with cr = yfhjJnuJj. This method is useful for solving 
problems similar to <[3j with nonlinear terms. 

XMDS2 permits the use of different spectral methods in 
each dimension. Figure [T] is an example of a lattice us- 
ing a Hermite-Gauss decomposition in the x dimension and a 
Fourier decomposition in the y dimension. As discussed in Sec- 
tion 3.4 the grid spacing is determined by the choice of spectral 



method. Full documentation of the spectral methods supported 
by XMDS2 and their uses is available from the XMDS2 website 

m. 

3.2. Fast spatial-to-spectral transforms 

In any nonlinear simulation, both the spatial and spectral rep- 
resentations of the solution will be required, as the problem will 
not be diagonal in either representation. Typically, the spa- 
tial representation is used for calculating the nonlinear terms, 
and the spectral representation for calculating derivatives. The 
two are linked by a linear transformation which can in general 
be performed with a matrix multiplication. The computational 
complexity of this operation is 0(N 2 ) for a single dimension. 
In higher dimensions, the use of a tensor product lattice (see 
Figure [T| enables the matrix multiplication to be factorised for 
each dimension. In two dimensions for example, the computa- 
tional complexity of a general spatial-to-spectral transformation 
is 0(N 2 N2 + NiN^). Without the use of a tensor product lattice, 
this cost would be 0(N 2 Nj). 

There are two cases in which we can reduce this computa- 
tional cost: when we can use the Fast Fourier Transform (FFT) 
algorithm, or when the basis functions alternate in parity. 

Spectral methods using complex exponentials, cosines or 
sines enable the use of the FFT algorithm and its variants 
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for transformations between spatial and spectral representa- 
tions. These cost only 0(N log N) in one dimension or 
0{N2N\ \ogN\ + N1N2 log N2) in two dimensions. 

If the basis functions have explicit, alternating parity 
X n (-x) - (-l)"X n {x) like the Hermite-Gauss functions, the Par- 
ity Matrix Multiplication Transform (PMMT) |[T4l can be used 
which is faster than a direct matrix multiplication in each di- 
mension. The idea is to separately transform the even and odd 
components of the solution, each of which costs O ((N/2) 2> j giv- 
ing a total cost of 0(N 2 /2). This factor of two reduction does 
not improve the overall scaling but can be a significant improve- 
ment for simulations dominated by the cost of the spatial-to- 
spectral transforms. 

3.3. Distributed memory parallelism 

The use of a tensor product lattice permits the problem to be 
parallelised by distributing a single dimension across the avail- 
able processes (see Figure |2j. The advantage of this method is 
that as the spatial-to-spectral transform can be factorised across 
different dimensions, when the problem is decomposed across 
the x dimension (as in Figure |2| a)), the transform over the y di- 
mension can be performed as a purely local operation to each 
process. 

To perform the spatial-to-spectral transform over the x di- 
mension, the problem must instead be decomposed across in 
the y dimension (as in Figure [2jb)). As simulations typically 
require spatial-to-spectral transforms to be performed over all 
dimensions, a distributed transpose operation is used to link 
different problem decompositions (see Figure [2]). This enables 
transforms to be performed over any dimension in a distributed 
simulation. 

3.4. Gaussian quadrature 

Gaussian quadrature is an exponentially accurate method for 
integrating functions. The key idea is to approximate 



J f{x) dxx^ f( Xi ) wt, 



(4) 



(a) 



(b) ,, 



Process 1 Process 2 Process 3 Process 4 
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Distributed transpose 



T 



where x, are the interpolation points and w, are weight factors. 
Gaussian quadrature takes advantage of the fact that the inter- 
polation points Xj do not need to be equally spaced. This means 
the 2N degrees of freedom {x,, w,} can be chosen to exactly in- 
tegrate 2N functions f(x), while it would only be possible to 
exactly integrate N functions if the w, were the only degrees of 
freedom. Further details about Gaussian quadrature are avail- 
able from flU [161 • 



Figure 2: An example of problem parallelisation on a tensor product lattice. 
The problem is distributed across the (a) x or (b) y dimensions. These two 
problem decompositions are linked by a distributed transpose operation. 



3.5. Method-of-lines explicit temporal integration 

In method-of-lines integration, each grid point is considered 
to have its own ODE and the problem is integrated as a system 
of coupled ODEs. XMDS2 employs a range of explicit integra- 
tion schemes for deterministic and stochastic problems: 

• semi-implicit method (deterministic order 2, stochastic or- 
der 1) ED, 
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• fourth-order Runge-Kutta (deterministic order 4, stochas- 
tic order 1/2) [TJ §3.7(v)], 

• ninth-order Runge-Kutta (deterministic order 9, stochastic 
order 1/2) d, 

• adaptive fourth-fifth order Runge-Kutta (deterministic 
only), (20) and 

• adaptive eighth-ninth order Runge-Kutta (deterministic 
only) HI. 

XMDS2's fixed-step method-of-lines integration methods 
support integrating stochastic differential equations that de- 
pend on Wiener (Gaussian) or jump (counting) processes. 
These stochastic differential equations must be entered in 
Stratonovich, not Ito form GTl . 

Although the fourth-order Runge-Kutta and ninth-order 
Runge-Kutta algorithms have lower order stochastic conver- 
gence than the semi-implicit method, we find that they can be 
useful for problems where the noise terms are a perturbation on 
the 'deterministic' dynamics. 

XMDS2 can run multiple paths (possibly distributed across 
multiple processors) to compute moments of the stochastic pro- 
cess. XMDS2 can also test the effect on the strong conver- 
gence Ell of the discretisation error of the propagation dimen- 
sion. This requires sampling the same stochastic trajectory with 
timesteps of multiple sizes. 

3.6. Interaction picture method 

The method-of-lines integration schemes are supported by 
the interaction picture method 11231 |24| which exactly solves a 
linear part of the differential equation. 

The idea is very similar to the interaction picture in quantum 
mechanics. The differential equation is split into two parts: a 
linear, exactly solvable component, and the remaining possibly 
nonlinear components. The differential equation is then trans- 
formed to remove the exactly solvable component. 

For a PDE of the form 



For example, for the nonlinear Schrodinger equation, 



df 
dt 



Z[f]+g{x,yJ), 



(5) 



where X is a linear operator that doesn't depend on time, the 
differential equation is transformed by defining the new quan- 



tity / — e ^'f which evolves as 



^=e-*g(x,y,e*< 



/)• 



(6) 



The new quantity / essentially has the simple dynamics due to 
X. removed. 

The interaction picture method is advantageous when the lin- 
ear operator X. has a faster characteristic timescale than the re- 
mainder of the differential system, which means that the func- 
tion / varies more slowly in time than the original /. This 
means that by solving the faster component separately and ex- 
actly, larger time-steps may be used on the remaining part of 
the differential equation while achieving the same solution ac- 
curacy. 



h 2 d 2 il/ 

m^ I = -—-^ + v(xW + u\ 

2m ox L 



difj 



(7) 



the spatial derivative term can have a faster characteristic 
timescale than the remainder of the system for high spatial res- 
olutions, corresponding to high momentum components. In the 
Fourier basis, the spatial derivative term in Eq. (j7]i becomes 



2 /,2 



h 2 k 

2m 



(8) 



If we do not use the interaction picture, the maximum value of 
k x increases linearly with the number of grid points, and the 
time step used must decrease as At oc 1 /N 2 in order to be able 
to resolve the evolution of those terms. The interaction picture 
method alleviates this problem by solving the spatial deriva- 
tive term exactly. Using the interaction picture method to solve 
for the evolution of the spatial derivative term enables (|7} to 
be solved with a time-step which is independent of the spatial 
resolution. 

The effect of the interaction picture method can be seen by 
solving (|7]i with an adaptive temporal integration method and 
comparing the number of time steps needed to achieve a given 
accuracy to that needed when calculating the derivatives explic- 
itly (but still using a spectral method). The results in Figure [3] 
demonstrate that the number of steps needed to solve the PDE 
using the interaction picture method is essentially independent 
of the spatial resolution, while for the explicit method, the num- 
ber of steps needed increases quadratic ally. 

The computational cost of the interaction picture method is 
small if the linear operator X is local in either the spatial or 
spectral basis. In this case, the application of e ±J -' to the quan- 
tity / can be calculated by transforming / to the appropriate 
basis (spatial or spectral), performing a local multiplication, 
and transforming back to the original basis. For fixed time-step 
algorithms, calculating the exponential function at every time 
step can be avoided by essentially redefining / at each time 
step so that only the quantities e ±jCA ' are needed. 

To make best use of the interaction picture method, the basis 
functions should be chosen to make all derivative terms local in 
the spectral basis. This ensures optimal scaling of the compu- 
tational effort with spatial resolution. 

Although the interaction picture method can be used with any 
integration software by applying the transformation manually, 
XMDS2 makes its use particularly easy by allowing the differ- 
ential equation to be entered in a form equivalent to <|5j with 
XMDS2 automatically making the transformation to |6]). This 
also enables easy comparisons to be made between the interac- 
tion picture and explicit methods. 

4. Examples 

In order to show the syntax of an XMDS2 script, as well 
as to demonstrate the ease with which simulations can be ex- 
tended, we consider the behaviour of a Bose-Einstein conden- 
sate (BEC) in a harmonic magnetic trap. 
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Figure 3: Comparison of the scaling of the interaction picture and 'explicit' 
methods with grid resolution for computing the evolution due to spatial deriva- 
tive terms. Both methods were used to integrate the PDE {5} with a fixed ac- 
curacy using an adaptive integrator. As the resolution is increased, the num- 
ber of steps remains approximately steady for the interaction picture method, 
while it increases quadratically (due to the second order spatial derivatives in 
l|9}) for the 'explicit' method. As the computational cost of each time step 
for both methods increases with resolution as 0(N log N) due to the use of 
Fourier transforms, the overall running time for the interaction picture method 
scales as 0(N log N) compared to 0(N 3 log AO for the 'explicit' method. The 
XMDS2 scripts used can be found in exampl.es/cpc_ip_scaling.xiiids and 
example/cpc_ex_scaling.xmds in the XMDS2 distribution. 



4.1. Example 1: Nonlinear Schrodinger equation 
( examples/cpc_examplel . xmds) 
Under a semiclassical approximation, the dynamics of the 
BEC will be governed by the nonlinear Schrodinger equation 
with a harmonic trapping potential. In dimensionless units this 
equation is written 



difr 



1 d 2 <p 
~2 dx* 



+ -]^iff+U\ifr\ 2 iff, 



(9) 



where x = y/mco/hx, t = cot, m is the atomic mass, to is the 
trapping frequency and U is the nonlinear energy in units of hco. 
XMDS2 is capable of solving much more complicated (sets) of 
PDEs, but this serves as a illustrative example. 

Our initial condition will specify the wavefunction at t — 0, 
and we choose 



t//(x,0) = V/V.T 1/4 exp(- 



■x 2 /2) 



(10) 



which is the ground state solution to Eq. ([9]) in the absence of 
the nonlinearity, normalized to TV atoms in total. 

We initially solve in one dimension, using a fourth-fifth or- 
der adaptive Runge-Kutta algorithm, evolving the system for a 
time t = 2n/co (one trap period), sampling 50 times and out- 
putting the real and imaginary parts of the wavefunction in po- 
sition space at every grid point. An XMDS2 script to solve this 
problem is shown in Figure |4] 

When XMDS2 is run on this script, it produces an optimized 
binary nonlinear_SE which is run to carry out the simulation. 
The result is shown in Figure [6] 

Changing parameters such as the domain or number of grid 
points, the number of sample points, integration interval, output 



moments, algorithm and precision used and so on is simply a 
matter of changing the contents of an XML tag, then re-running 
XMDS2 on the script to produce the new executable. While it 
is trivial to change such parameters, it is also easy to extend the 
simulation in more complex ways. 

4.2. Example 2: Higher dimensions 
( examples/ cpc_example2. xmds) 

If one wished to run the simulation in two dimensions rather 
than one, all that is required is adding the element 

<dimension name="y" lattice="512" 

domain="(-7, 7)" /> 

to the transverse_dimensions element, changing the basis="x" 
attribute of the sampling_group element to basis="x y", adding a 
"0 . 5*y*y" term to the potential "V" and initial condition, and adding 
a "-i*0 . 5*ky*ky" term to the kinetic energy operator "T". 

4.3. Example 3: Different transforms 
( examples/ cpc_example3 . xmds) 

This problem is obviously symmetric about x = 0, so it is a waste 
of computational resources to simulate the problem on both sides of 
the origin. Since the differential equation and the boundary conditions 
are symmetric, by using the discrete cosine transform rather than the 
default exponential Fourier transform, we need only carry out the sim- 
ulation on half the interval, and use only half the number of grid points 
for the same accuracy. This is accomplished simply by changing the 
content of the transverse_dimensions element to be 



<dimension name="x' 
domain=' 



lattice="256" 
(0, 7)" transform 



"dct"/> 
with 



MPI 



4.4. Example 4: Easy parallelization 
( examples/ cpc_exampleJf.. xmds) 

As this is a deterministic simulation, if it has two or more di- 
mensions, one can parallelize the simulation simply by adding the 
<driver name="distributed-mpi" /> tag to the script. This 
would result in a binary that could be run across, for example, four 
CPUs with the command 

mpirun -n 4 nonlinear_SE 

The fact that two- or higher-dimensional deterministic simulations, 
as well as stochastic simulations of any dimension, can be trivially 
parallelized using a single line in a script, without spending days (or 
weeks) writing and debugging bespoke code, is one of XMDS2's most 
powerful features. The runtime scaling of this simulation with the 
number of processes is illustrated in Figure[5] 

4.5. Example 5: Non-local terms 
( examples/ cpc_example5 . xmds) 

Many problems will involve non-local interactions that occur in the 
form of a convolution j f(r - r')g(r') dr' . For example, within the 
context of the current problem, if the BEC were charged there would 
be an additional potential of the form 

„2 7 2 



V(X): 



e-Z- 
4tt<e 



J \x-x'\ 



)\ dx 



(11) 



where eZ is the charge associated with each particle. While this term 
could be explicitly integrated within XMDS2 it is more efficient to 
make use of convolutions and the speed of fast Fourier transforms. 
This is done using <computed_vector> elements, which are de- 
scribed in detail on our website |T). 
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simulation xmds-version="2"> 
<name> nonlinear_SE </name> 



<f eatures> 
<globals> 
<! [CDATA[ 

real N = 10.0; // number of atoms 
real g = 1.0; // nonlinear coupling 

]]> 
</globals> 
</f eatures> 

<geometry> 

<propagation_dimension> t </propagation_dimension> 
<transverse_dimensions> 

<dimension name="x" lattice="51Z" 
domain="(-7, 7)" /> 
</transverse_dimensions> 
</geometry> 



<vector name="potential" type="real"> 
<components> V </components> 
<initialisation> 
<! [CDATA[ 

V = 0.5*x*x; 
]]> 

</initialisation> 
</vector> 



<vector name="wavef unction" type="complex"> 

<components> psi </components> 
<initialisation> 
<! [CDATA[ 

psi = sqrt(N) * pow(M_PI, -0.Z5) * exp(-x*x/2); 
]]> 

</initialisation> 
</vector> 



Simulation constants 



Define the geometry 

The potential is real 
Define the potential: 



V(x) = -x 2 



The wavefunction is complex 

Define the wavefunction initial condition: 

ip(x, t = 0) = VNn- 1/4 exp(-x 2 /2) 



<sequence> 

<integrate algorithm="ARK45" interval="6.28" 

tolerance="le-5"> 
<samples> 50 </samples> -4 



<operators> 

<operator kind="ip" > 

<operator_names> T </operator_names> 
<![CDATA[ 

T = -i * 0.5 * kx * kx; 

]]> 
</operator> 

<integration_vectors> 

wavefunction 
</integration_vectors> 
<dependencies> potential </dependencies> 
<! [CDATA[ 

dpsi_dt = T[psi] - i * V * psi 

- i * g * modZ(psi) * psi; 

]]> 

</operators> 
</integrate> 
</sequence> 



<output format="hdf 5"> ' 

<sampling_group initial_sample="yes" basis="x"> 
<dependencies> wavefunction </dependencies> 
<moments> psi real psiimag </moments> 
<![CDATA[ 

psi real = Re(psi); 
psiimag = Im(psi); 
]]> 

</ sampl i ng_group> 
</output> 
</simulation> 




Integration method, interval and tolerance 

Sample output 50 times 

Define the interaction picture operator: 

Evolve the wavefunction in time 
We need the potential below: 

Define the evolution equation: 

— = T[ip] - iV(x)iji - ig |V>| 2 i> 

Output in HDF5 format 
Sample at f=0 

Sample the spatial representation 



Output the real and imaginary 
parts of tp 



Figure 4: Annotated example XMDS2 script for integrating equation |5J. This script can be found in examples/cpc_examplel . xmds in the XMDS2 distribution. 
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.6. Example 6: 

( examp I es/cp c_examp I e6 . xmds) 



Stochastics 



Cl 

■a 
a) 
<x> 

CL 

CO 



12 
10 

8 
6 
4 
2 



Example 4 (2D) performance 
-Perfect scaling 




6 8 
MPI Processes 



70 



(b) 



Example 4 (3D) performance 
Perfect scaling 




40 

MPI Processes 



80 



Figure 5: Example runtime scaling using MPI on (a) a single com- 
puter, and (b) a supercomputer. Figure (a) demonstrates the simulation 
examples/cpc_example4 . xmds run on a Linux computer with two Xeon 
5675 CPUs running at 3.07GHz. Each CPU has 6 execution cores. Figure 
(b) demonstrates the simulation examples/cpc_example4_3D . xmds run on 
the NCI National Facility supercomputer, 'vayu'. The modified simulation is 
extended to three dimensions with 256 points in each to demonstrate perfor- 
mance on larger problems. Note that optimal parallelisation for these problems 
is achieved when the number of grid points in the first dimension (256) is divis- 
ible by the number of processes. 



As a final tweak to this example, we will make use of XMDS2's 
stochastic features to add noise. XMDS2 has a number of fast random 
number generators built in, which are capable of producing Gaussian, 
Poissonian and uniform probability distributions, and applying them as 
Wiener or jump processes during stochastic integration. This enables 
the simulation of stochastic differential equations, which are useful in 
fields such quantum field theory, mathematical finance, and many oth- 
ers. For this example we will simply use noise to model perturbations 
of the magnetic trap — that is, the trapping potential will be slightly 
noisy, due to it moving around. To do this we define a noise vector 

<noise_vector name="trapNoise" kind="wiener" 
type="real" method="dsfmt"> 
<components> noise_x </components> 
</noise_vector> 

change the potential term in the equation of motion in the script to 

- i * (V + g * mod2(psi) + alpha*noise_x) * psi 

where alpha is a constant governing the magnitude of the noise, 
and add a <dependencies>trapNoise</dependencies> tag to the 
<initialisation> block of the potential vector. This would add 
a time-dependent Gaussian-distributed noise to the potential. If we 
wished to average over many different realisations of this noise, we 
could add the tag 

<driver name="mpi-mrilti-path" paths="100" /> 

to the script, which would run the simulation 100 times, and average 
over whichever results were requested in < output > section. Such a 
simulation could be trivially run over any number of CPU cores with 
near perfect scaling. 

Note that the mpi-multi-path driver should only be used for 
stochastic simulations where individual realisations are independent, 
in contrast, the distributed-mpi driver parallelises a single de- 
terministic simulation. As the different components of a determin- 
istic simulation will in general be coupled, the distributed-mpi 
driver necessarily incurs a larger communication overhead than the 
mpi-multi-path driver. In general, the distributed-mpi driver 
can be used to parallelise a single realisation of a stochastic simula- 
tion, but if many paths are needed, the mpi-multi-path driver will 
be preferable. 



5. Software used 

XMDS2 makes use of the following external libraries in its gener- 
ated C++ simulations: 

• FFTW3 1 25 1 for FFTs and MPI distributed transpose operations, 

• dSFMT 1261 for random number generation, 

• MPI for inter-process communication, 

• HDF5 for data input and output, and 

• GNU Scientific Library for special function evaluation. 

XMDS2 itself also uses the following Python libraries when generat- 
ing simulations: Cheetah, pyparsing, lxml, h5py, mpmath, and numpy. 
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Figure 6: Solution to the nonlinear Schrodinger equation given by Eq. |9j. The 
density \i//(x, r)| 2 is shown evolving over one trap period. 

6. Conclusion 

Using XMDS2 for simulations accelerates development time, pro- 
duces code that executes extremely quickly, and also produces a self- 
documenting workflow, as output data is wrapped with the compact 
XML code used to produce it. 

XMDS2 particularly excels at providing a smooth transition from 
a low-dimensional simulation to a higher-dimensional one, from a de- 
terministic simulation to a stochastic one, or from a single-processor 
simulation to a distributed simulation running in parallel across mul- 
tiple computers (or on a supercomputer). In hand-written codes, un- 
less they were initially written with such a potential future extension 
in mind, each such change would require significant effort in rewrit- 
ing the code. In XMDS2 such changes require only minimal change 
to the input script. This encourages users to create test simulations 
of a simpler system (e.g. reduced dimensionality), which makes the 
code run faster, allowing problems in the input script to be found and 
fixed more quickly. Later, the simulation can be scaled up to the full 
problem. Fundamentally, the ease with which codes can be generated 
encourages experimentation with different types of simulations, as the 
time taken to create the code is no longer the rate-limiting factor. 

The installers, documentation and examples for XMDS2 can be 
found at the website [ 1 ]. This same documentation is available in the 
documentation/ directory of the XMDS2 distribution. 
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